# TODO: create PM image from aod and temp (PM regression)
# 
# Author: phamha
###############################################################################

#import library
library(gstat)
library(RPostgreSQL)
library(stringr)
library(raster)

#import source
setwd("/home/phamha/project/Regression")
source("constant.R")
source("util.R")

data_folder = "/var/www/html/"
met_folder = "/home/phamha/MetTiff/"
myd_folder = "/var/www/html/fimo/apom/Product/MYD01PM/"
mod_folder = "/var/www/html/fimo/apom/Product/MOD01PM/"

mod_time_query="select aqstime from apom.satresampmod04 where filename like '%"
myd_time_query="select aqstime from apom.satresampmyd04 where filename like '%"
mod_query = "SELECT mod04.aqstime, mod04.filename,mod04.filepath,mod07.filename as temp_filename,mod07.filepath as temp_filepath FROM apom.satresampmod04 as mod04 inner join apom.satresampmod07temperature as mod07 ON (mod04.aqstime = mod07.aqstime) where mod04.aqstime < '2014-09-12 00:00:00'::timestamp"
myd_query = "SELECT mod04.aqstime, mod04.filename,mod04.filepath,mod07.filename as temp_filename,mod07.filepath as temp_filepath FROM apom.satresampmyd04 as mod04 inner join apom.satresampmyd07temperature as mod07 ON (mod04.aqstime = mod07.aqstime) where mod04.aqstime < '2014-09-12 00:00:00'::timestamp"

create_regress_pm = function(sate_data,model_type,cook_type,aod_file,temp_file){
	if(sate_data=="mod"){
		time_query = mod_time_query
		min_aod = min_aod_mod
		max_aod = max_aod_mod
		min_temp = min_temp_mod
		max_temp = max_temp_mod
		
	}else{
		time_query = myd_time_query
		min_aod = min_aod_myd
		max_aod = max_aod_myd
		min_temp = min_temp_myd
		max_temp = max_temp_myd
	}
	
	start_index = regexpr("M[^M]*$",aod_file)
	end_index = regexpr(".tif",aod_file) - 1
	file_name = substr(aod_file,start_index,end_index)
	
	time_query = paste(time_query,file_name,"%'",sep="")
	data = getDataFromDB(time_query)
	
	mod04_aqstime = data$aqstime[1]
	aqstime = strptime(mod04_aqstime,format="%Y-%m-%d %H:%M:%S")
	aqstime = aqstime + 25200	
	month = format.Date(aqstime,"%m")
	year = format.Date(aqstime,"%Y")

	aod_dataset = raster(aod_file)
	aod_data = values(aod_dataset)
	aod_data = aod_data * 0.00100000004749745
		
	corxy = coordinates(aod_dataset)
	x = corxy[,'x']
	y = corxy[,'y']
		

	temp_dataset = raster(temp_file)
	temp_data = values(temp_dataset)
	temp_data = (temp_data + 15000) * 0.00999999977648258
		
	avg_temp_file  = paste(met_folder,"temp",as.numeric(month),".tif",sep="")
	avg_rh_file    = paste(met_folder,"rh",as.numeric(month),".tif",sep="")
	avg_preci_file = paste(met_folder,"preci",as.numeric(month),".tif",sep="")
		
		
	avg_temp_dataset  =  raster(avg_temp_file)
	avg_rh_dataset    =  raster(avg_rh_file)
	avg_preci_dataset =  raster(avg_preci_file)
		
	avg_temp_data  = values(avg_temp_dataset) 
	avg_rh_data    = values(avg_rh_dataset) 
	avg_preci_data = values(avg_preci_dataset)
		

	pm25 = regress_predict(sate_data,model_type,cook_type,aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data)

	table = data.frame(x,y,aod_data,temp_data,avg_temp_data,avg_rh_data,avg_preci_data,pm25)
	table$pm25[table$aod_data<min_aod|table$aod_data>max_aod|table$temp_data<min_temp|table$temp_data>max_temp]<-NA
	og_raster = aod_dataset
	totalCell = ncell(og_raster)
	og_raster[1:totalCell] = table$pm25

	pm_file = str_replace(aod_file,".tif","_.tif")
	writeRaster(og_raster,filename=pm_file,format="GTiff")
	print(file_name)
			
}
aod_file = "C:/MOD04_L2.A2014024.0310.051.2014028020029_rg.tif"
temp_file = "C:/MOD04_L2.A2014024.0310.051.2014028020029_rg.tif"
create_regress_pm("mod","linear","4np",aod_file,temp_file)
print("done")
